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ABSTIfACT 

We present a consensus-based distributed particle filter (PF) for 
wireless sensor networks. Each sensor runs a local PF to compute a 
global state estimate that takes into account the measurements of all 
sensors. The local PFs use the joint (all-sensors) likelihood function, 
which is calculated in a distributed way by a novel generalization of 
the likelihood consensus scheme. A performance improvement (or a 
reduction of the required number of particles) is achieved by a novel 
distributed, consensus-based method for adapting the proposal den- 
sities of the local PFs. The performance of the proposed distributed 
PF is demonstrated for a target tracking problem. 

Index Terms — Distributed particle filter, likelihood consensus, 
distributed proposal density adaptation, target tracking, wireless sen- 
sor network. 

1. INTRODUCTION 

We consider distributed sequential state estimation in a wireless sen- 
sor network. For general nonlinear/non-Gaussian scenarios, the par- 
ticle filter (PF) is often the estimation method of choice 1 1 1. In this 
paper, extending our work in | 2~4 |, we propose a distributed PF that 
uses a novel distributed scheme for proposal density adaptation. As 
in (2fi4|, each sensor runs a local PF that computes a global state 
estimate incorporating the measurements of all sensors. The local 
PFs use the joint (all-sensors) likelihood function (JLF), which is 
computed in a decentralized way by means of the likelihood con- 
sensus (LC) scheme. Here, we present a generalized form of the 
LC originally proposed in [2|, which is suited to a general measure- 
ment model (i.e., it is not limited to additive Gaussian measurement 
noises l2ll3l or likelihoods from the exponential family [4|). 

Our main contribution is a novel distributed, consensus-based 
scheme for adapting the proposal densities (PDs) used by the local 
PFs. Adapted PDs can yield a significant performance improvement 
or, alternatively, a significant reduction of the required number of 
particles |5 |. In our adaptation scheme, local PDs computed by the 
individual sensors are fused in a distributed way by means of con- 
sensus algorithms, thereby providing to each local PF a global PD 
reflecting all measurements. To make our scheme computationally 
feasible, we use Gaussian approximations for the local and global 
PDs. Our PD adaptation scheme differs from that proposed in O in 
that it is distributed and it uses Gaussian PD approximations. 

Consensus-based distributed PFs with PD adaptation have been 
recently proposed in |7 8|. The distributed PD adaptation scheme 
of 1 7 1 employs min- and max-consensus to construct a set capturing 
most of the posterior probability mass. This set is used to calcu- 
late a distorted state-transition density, which serves as PD. Our dis- 
tributed PD adaptation scheme has a lower complexity than that of 
^|. The communication requirements of our PD adaptation scheme 
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are somewhat higher, but the overall communication requirements of 
our distributed PF are still much lower than those of the distributed 
PF of |7|, and simulation results demonstrate a better estimation per- 
formance of our distributed PF. In | 8|, a distributed unscented PF is 
presented. This method employs a PD adaptation which, however, 
is not distributed: the PD used at each sensor is only based on the 
local measurement. Again, simulation results demonstrate a better 
estimation performance of our distributed PF, which however comes 
at the cost of higher communication requirements. 

This paper is organized as follows. In Section [2l we introduce 
the system model and review the principles of sequential Bayesian 
estimation. The LC-based distributed PF and the new generalized 
LC scheme are described in Section [3] In Section |4l we present 
the proposed distributed PD adaptation scheme. Finally, Section [5] 
reports simulation results for a target tracking problem. 

2. SEQUENTIAL BAYESIAN STATE ESTIMATION 

We consider a random, time-varying state vector x„ = {xn,i ■ ■ ■ 
Xn^Ai)^- The state evolves according to the state-transition model 



gn(x„_l, Un 



1,2, 



(1) 



where u„ is white driving noise with a known probability density 
function (pdf) /(u„). At time n, x„ is sensed by a sensor network 
with K sensors according to the measurement models 



i(x„, V„ 



= 1,2, 



(2) 



Here, z„,fe of dimension Nn^k is the measurement at time n and at 
sensor k, and 'Vn,k is measurement noise with a known pdf /(v„,fc). 
We assume that (i) v„,fc and v„/ are independent unless (n, k) — 
{n', k'); (ii) the initial state xo and the sequences u„ and v„,fc are 
all independent; and (iii) sensor k knows gn(-, •) and h„,fe(-, ■) for 
all n, but it does not know h„.j;/ (•, •) for k'^ k. 

The state-transition and measurement models (TJ and ^ together 
with our statistical assumptions determine the state-transition pdf 
/(x„|x,i_i), the local likelihood function /(z,j.fc |x,i), and the JLF 
/(z,i|x„). Here, z„ = {zZ,i ' ■ ' zj /.^)^ denotes the vector contain- 
ing all sensor measurements at time n. Due to ^ and the indepen- 
dence of all v„^fc, the JLF is given by 



/(z„|x„) = ]^/(z„,fc|x„) 



(3) 



Our goal is estimation of the state x„ from all sensor measure- 
ments from time 1 to time n, zi;„ = (z^- ■ • zX)^- To this end, we 
consider the minimum mean-square error (MMSE) estimator |9| 



E{x„|zi;„} = / x„/(x„|zi;„) dx„ 



(4) 



The posterior pdf /(x,i|zi:„) in Q can be calculated sequentially 
from the previous posterior /(x„_i |zi;„_i) and the JLF /(z„ix„) 
\Wi . A computationally feasible approximation to this sequential 



MMSE state estimation is provided by the PF, which represents the 
posterior pdf /(x„ |zi;„) by a set of weighted particles 1 1 1. 

3. LC-BASED DISTRIBUTED PARTICLE FILTER 

The proposed LC-based distributed PF (LC-DPF) differs from our 
previous work in |2-4| by the generalized LC and the distributed 
PD adaptation (presented in Sections 13.21 and |4] respectively). As 
in 0^1, each sensor tracks a particle representation of the global 
posterior /(x„jzi:,i) using a local PF. At each time n, each local 
PF obtains a state estimate x„,fc that is based on zi-.n, i.e., all sen- 
sor measurements up to time n. This requires knowledge of the JLF 
/(z„ |xn) as a function of x„. An approximation of the JLF is pro- 
vided to all sensors in a distributed way by means of the generalized 
LC. No communication between distant sensors or complex rout- 
ing protocols are required. Also, no particles, local state estimates, 
or measurements are communicated between the sensors. The pro- 
posed distributed PD adaptation scheme can yield a significant per- 
formance improvement or, alternatively, a significant reduction of 
the number of particles (and, thus, of the computational complexity); 
this comes at the cost of an increase in inter-sensor communications. 



Initialization: The recursive procedure defined by Steps 1-7 is 
initialized at time n = by J particles Xq"*^ randomly drawn from an 
appropriate prior pdf /(xq), and by equal weights WqI = 1/J. 

3.2. Generalized Likelihood Consensus 

We now present the generalized LC scheme that is used in Step 5 
to provide an approximate JLF to each sensor. In contrast to our 
previous work |2-4 |, this scheme is not limited to likelihoods with 
exponential form or to additive Gaussian measurement noises; it is 
suitable for any type of likelihood and any measurement model l|2j. 
To derive the generalized LC, we first take the logarithm of l|3j: 



log/(z„|x„) = ^log/(z„,fc|x„) 



(6) 



Unfortunately, a consensus-based distributed calculation of l[6j is not 
possible in general because the terms of the sum depend on the un- 
known state Xn. We therefore use the following approximate (finite- 
dimensional) basis expansions of the local log-likelihoods: 



3.1. Local PF Algorithm 

At a given time n > 1, the local PF at sensor k performs the following 
steps, which are identical for all k: 

Step 1: A resampling (1 1 is applied to the J particles {x^^^i fc}j_i 
with corresponding weights k]'' -\ (calculated at time n— 1) 
that represent the previous global posterior /(x„_i |zi:„_i) at sen- 
sor k. This produces J resampled particles {x^^^i k]^j-\- 

Step 2: Temporary particles {x^^^^}'^_j^ are randomly drawn 



from /(x„|x. 



-l.fc) = /(Xn|x„_l) 

sian approximation J\f{'x.n', fJ^'n^ky CJ, j. 



/(x„|zi:„_i) is calculated according to 
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of the "predicted posterior" 
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J = l 



(5) 

Step 3 (jointly performed by all sensors, using communication 
with neighboring sensors): An adapted Gaussian PD g(x„; z„) = 
A/'(xn; /.tn, Cn) iuvolviug all sensor measurements is computed 
from the /xj^ C'^ f.,, and z„ k' (k' = 1, • • • , K) by means of the 
distributed, consensus-based scheme described in Section|4] 

Step 4: J particles {x^^].} are drawn from the PD ^(xn; z„). 

Step 5 (jointly performed by all sensors, using communication 
with neighboring sensors): An approximation /(z„|x„) of the JLF 
/(z„|x„) is computed in a distributed way by means of the gener- 
alized LC described in Section [J!2l using the particles 
drawn in Step 4. 

(7) 

Step 6: Weights associated with the particles x^ j. drawn in Step 
4 are calculated according to 



„0) 



U) |„(i) 



.,J, 



where 7 is chosen such that 



„(i) 



U) 



}^,_^, an approximation of the global 



Step 7: From {xj^-' 
MMSE state estimate x^^^^ in l|4j is computed according to 

J 

E(i) (i) 
'^n,k '^n.k ■ 



log./(z 



(7) 



Here, a„,fc.r(z„,fc) are expansion coefficients that contain all sensor- 
local information (including the sensor measurement z„_fc) and 
<(5n,i (x„) are fixed, sensor-independent basis functions that are 
assumed to be known to all sensors. Substituting into (|6}, we 
obtain 

R 

log/(z„|x„) ^ ^ an,r(z„) V9„,r(x„) , (8) 
r = l 

With 

K 

a7i,r(z„) = ^ an,fe,r(z„,fe) . (9) 
k = l 

The sum over all sensors in ^ can be easily computed in a dis- 
tributed way by means of a consensus algorithm 1 1 Ij since the terms 
of the sum are real numbers (not functions of x,i). 

By exponentiating ([8), we finally obtain the following approxi- 
mation of the JLF, denoted /(z„|x„): 

/(z„[x„) « /(z„|x„) = exp^^a„,r(z„)v'„,r(x„)^ . (10) 

Therefore, a sensor that knows the coefficients a„.r(z„) is able to 
evaluate the approximate JLF /(z„ |x„) for all values of x„. In fact, 

the vector of all coefficients, a„(z„) = (a„,i(zn) ■ ■ ■ a„,R{z„)Y, 
can be viewed as a sufficient statistic 1 9 1 that epitomizes the total 
measurement Zn within the limits of the approximation The ex- 
pressions ^ and l llOt allow a distributed, consensus-based calcula- 
tion of /(z„|x„) due to the following key facts, (i) The coefficients 
a.n,r{z„) do not depend on the state x„ but contain the information 
of all sensors (i.e., the expansion coefficients an,fe,r(z„,fe) for all fe). 

(ii) The state x„ enters into /(z„|x„) only via the basis functions 
'Pn,r{-), which are sensor-independent and known to each sensor. 

(iii) According to the a„,r(zn) are sums in which each term 
contains only local information of a single sensor. 

At each time n, the expansion coefficients cir,n.fe(z„,fc) in ^ are 
calculated locally at each sensor k by means of least squares fit- 



ting (T2I based on the J data points {x^^j, , log f{z„,k |x„ j.) } 



Here, the use of the particles x^-''^, drawn in Step 4 of the local PF 
algorithm ensures a good approximation in those regions of the state- 



space where the approximate JLF /(z„|x„) is evaluated in Step 6. 
(We assume that /(z„,fc |x^^^^j^) / 0.) 

The steps of the generalized LC scheme performed at a given time 
n can now be summarized as follows: 

Step 1: Sensor k calculates the expansion coefficients Q:r,n,fc(z,i.fc) 
in ([7} using least squares fitting. 

Step 2: The coefficients Qr.n.fc (zn.fe) of all sensors k are added in 
a distributed way using a consensus algorithm. One instance of that 
algorithm is employed for each r G {1, . . . , R}; all instances are ex- 
ecuted in parallel. After a sufficient number of consensus iterations, 
the ar,n(zn) (see ^) for all r are available at each sensor. 

Step 3: Using the ar.n{zn), each sensor is able to evaluate the 
approximate JLF /(z„ |x„) for any value of x„ according to JlOt . 

The proposed LC-DPF (without the PD adaptation described in 
Section |4l( requires the transmission of IR real numbers by each 
sensor at each time n, where / is the number of consensus iterations 
performed by each consensus algorithm and 7? (cf . ([7}) is the number 
of consensus algorithms executed in parallel. All transmissions are 
to neighboring sensors only, and their number does not depend on the 
measurement dimensions Nn^k- Thus, the LC-DPF is particularly 
attractive in the case of high-dimensional measurements. 

4. DISTRIBUTED PROPOSAL ADAPTATION 

We now present our distributed scheme for calculating the adapted 
PD g(x,i; z,i) (Step 3 of the local PF algorithm in Section [3Tb . This 
scheme can be summarized as follows. First, a "pre-distorted" local 
posterior is calculated at each sensor. The local posteriors are then 
fused via a distributed fusion rule to obtain a global posterior, which 
is used as PD by each local PFQ This PD takes into account the 
measurements of all sensors, which is appropriate in view of the fact 
that the JLF is used in Step 6. Our approach is inspired by the one 
from 1 13], which however was proposed in a different context and 
uses a fusion rule different from ours. 

We first note that the global posterior /(x„|zi;„) can be written 
(up to a normalization factor) as 

/(Xn|zi;„) = /(x„|zi;„_l,Z„) 

OC /(z n |Xn , Zl;n — 1 ) f ip^n |zi;n — 1 ) 
= /(Z,i|x„) /(X„|zi;„_l) 
■ K 

Y[ /(z„,fc|x„ 



/(X„|zi:„_l) . (11) 



Let us suppose that each sensor k calculates a (pre-distorted, non- 
normalized) local pseudoposterior defined as 

/(X„|Z1:„_1, Z,i,fc) = /(Z„,fc|x„) /^^"'^(X„|zi:,i_l) . (12) 

The product of all local pseudoposteriors equals the global posterior 
up to a factor: 



n/(> 



Zl:, 



-l,Zn,fc) ~ 



K 



n /(^" 



OC /(X„|Z1;„) , 



/(X„|zi:„_l) (13) 



where Jilt has been used. This posterior reflects all sensor mea- 
surements and could be employed as the global PD. However, for 
a simple distributed computation of ( |13t , we use Gaussian approx- 
imations of the local pseudoposteriors and the global posterior, 
i.e., /(x„|zi:„_i, z„,fc) f» A/'(x„; /i„,fe, C„,fc) and /(x„|zi;„) f» 



'Note that the global posterior used as PD is different from the global 
posterior that is obtained by the PF as described in Section[3] 



g(x„;z,i) — A/'(x„; /!„, C„). Then, using ( IT3t and the rules 
for a product of Gaussian densities |14), we obtain the following 
expressions of the mean and covariance of the PD g(x„; z,i): 



— Cn C^ j. p,„,k 



5Z 



(14) 



The sums over all sensors in these expressions can be easily calcu- 
lated in a distributed way using consensus algorithms. 

To calculate the Gaussian approximation A/'(x„; fin.k, C„,k) of 
the local pseudoposterior /(x„|zi:,i_i, z„,fe), we note that l ll2b 
is the measurement update step of a Bayesian filter using the 
pre-distorted predicted posterior /^''^'^(x„|zi:„_i) instead of the 
true predicted posterior /(x„ |zi;„_i). Furthermore, each sen- 
sor calculated a Gaussian approximation of the predicted poste- 
rior, /(x„|zi;„_i) f» A/'(x„;/i^_fe,C^_fe) (see Step 2 in Section 

13. Il l: this entails the Gaussian approximation /^''^(x,i|zi:„_i) « 
A/'(x„; A'Cjj^). Since Gaussian models are thus used for 

both /(x„|zi;„-i, z„,fc) and /^''^x„|zi;„-i), we propose to per- 
form the measurement update in l ll2b by means of the update step of 
a Gaussian filter il0lll5Hl71 . This is done locally at each sensor. 

The operations of the proposed PD adaptation scheme performed 
at time n can now be summarized as follows: 

Step 1: Each sensor k computes the mean p,n,k and covariance 
C„,fe of the Gaussian approximation of the local pseudoposterior, 

Af{-x.n; fin,k,Cn,k) ~ / (x„ |zi 1 , z„,fc ) . Thls Is douc locally 
by performing a Gaussian filter update step with input mean /xlj.^, 
input covariance KC'^j^, and measurement z„,fc. Here, /i^.fe and 
C'„ f. were obtained locally according to (|5}. 

Step 2: Consensus algorithms are used to calculate the sums over 
all sensors in J14t . This step requires communication with neigh- 
boring sensors. In total, / [M + M{M -1- 1)/2 -|- 1] real numbers are 
communicated by each sensor at time n. Here, I denotes the number 
of consensus iterations and AI denotes the dimension of the state x„ . 
After convergence of the consensus algorithms, each sensor obtained 
the global PD <j(x„; z„) = A/'(x,i; C„). 

5. SIMULATION RESULTS 

We consider a target tracking application using acoustic ampli- 
tude sensors. The target is represented by the vector t„ — 
(xn Un in yn)^ Containing the target's 2D position and 2D velocity 
in the x-y plane. The vector t,i evolves with time n according to 
T„ = Gt„_i +Wu^, n = 1, 2, . . . , where the matrices G G R**^* 
and W G M}^'^ are chosen as in |4| and the are independent 
and identically distributed according to ^ A/'(0, C^/) with 
C^/ — diag(0. 0033, 0.0033). The target motion model specified 
above is however assumed unknown to the simulated PFs. There- 
fore, all simulated PFs use a random walk model x„ = x„_i + u„, 
where the state x„ = {xn y-n)^ represents the position of the target 
and u„ ^Af{0, C„) with C„ = diag(0.0528, 0.0528) (cf. GJ). 

The target emits a sound of constant amplitude A — 10, which 
is sensed by acoustic amplitude sensors. The (scalar) measurement 
z„,k of sensor k is given by (cf. (|2}) 



A 



Zn,k = 



+ V„,k 



where ^^.k is the position of sensor k at time n and Vn,k ~ A/'(0, a^) 
with — 0.00005. This value of cr^ yields a peaky likelihood, 
which highlights the performance gains of PD adaptation. The net- 
work consists of K = 25 sensors that are deployed on a jittered grid 
within a rectangular region of size 40 m x 40 m. Each sensor com- 
municates with other sensors within a range of 18 m. 




For LC, unless stated otherwise, we approximate log /(z„,fc |x„) 
by a multivariate polynomial of degree Rp = 6; this leads to a basis 
expansion (O of degree R = {^"j^") = 28. The sums in ^ and (O 
are computed by / = 15 iterations of an average consensus algorithm 
with Metropolis weights 1 18|. For PD adaptation (Step 1 in Section 
0, the update step of an unscented Kalman filter 1 15| is used. 

We compare the proposed LC-DPF with the distributed PFs pre- 
sented in |7| and |8| (referred to as DPF-1 and DPF-2, respectively) 
and with a centralized PF (CPF) that processes all sensor measure- 
ments at a fusion center. The CPF uses an adapted PD that is com- 
puted using a (centralized) unscented Kalman filter. The number 
of particles at each sensor of the distributed PFs and at the fusion 
center of the CPF is J = 200 unless stated otherwise. As a per- 
formance measure, we use the root-mean-square error of the state 
estimates Xn.fc, denoted RMSE„, which is computed as the square 
root of the average of the squared estimation error over all sensors 
and over 1000 simulation runs. We also compute the average RMSE 
(ARMSE) by averaging RMSE^ over all 200 simulated time instants 
n and taking the square root of the result. 

Fig. [T] shows the temporal evolution of RMSEn. It can be seen 
that the performance of LC-DPF is almost as good as that of CPF 
and better than that of DPF-1 and DPF-2. The communication re- 
quirements of LC-DPF are lower than those of DPF-1 but higher 
than those of DPF-2: the total counts of real numbers transmitted 
by LC-DPF, DPF-1, and DPF-2 during one time step in the entire 
network (all sensors) are 12375, 76875, and 1875, respectively. 

Fig.|2]shows the ARMSE versus the degree Rp of the polynomial 
used to approximate the local log-likelihood functions. As expected, 
the ARMSE of LC-DPF decreases with growing Rp and approaches 
that of CPF. Note, however, that the communication requirements 
increase with growing Rp. 

Finally, Fig.[3]shows the dependence of the ARMSE on the num- 
ber J of particles for LC-DPF and an LC-DPF without proposal 
adaptation (abbreviated as LC-DPF-NA) |4|. As we can see, the 
ARMSE of LC-DPF is significantly lower than that of LC-DPF-NA, 
even if J is large. This demonstrates the performance improvement 
achieved by our distributed PD adaptation scheme. 

6. CONCLUSION 

We presented a consensus-based distributed particle filter (PF) for 
wireless sensor networks. The state estimates computed by the local 
PFs at the various sensors reflect the past and present measurements 
of all sensors. This is enabled by a generalized likelihood consen- 
sus scheme, which performs a distributed approximate calculation of 
the joint likelihood function for general measurement models. Our 
main contribution was a distributed method for adapting the proposal 
density (PD) used by the local PFs. This method is based on a Gaus- 
sian model for the PD, whose mean and covariance are computed 
in a distributed way by consensus algorithms and by the update step 
of a Gaussian filter. Simulation results demonstrated the good per- 



formance of the proposed distributed PF and the large performance 
gains achieved by the proposed PD adaptation method. 
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